// BFSTEST : Test breadth-first search in a graph.
// 
// example: cat sample.txt | ./bfstest 1
//
// John R. Gilbert, 17 Feb 2011

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <signal.h>
#include <sys/time.h>
#include <omp.h>

#define min(a,b) ((a) <= (b) ? a : b)
#define CHUNKSIZE 100

#define INTERVAL 500
#define US_INTERVAL 0 

#include <unistd.h>
#include <stdarg.h>

#define PNG_DEBUG 3
#include <png.h>

typedef struct graphstruct { // A graph in compressed-adjacency-list (CSR) form
    int nv;            // number of vertices
    int ne;            // number of edges
    int *nbr;          // array of neighbors of all vertices
    int *firstnbr;     // index in nbr[] of first neighbor of each vtx
} graph;


void write_png_file(char* );

void abort_(const char * s, ...)
{
        va_list args;
        va_start(args, s);
        vfprintf(stderr, s, args);
        fprintf(stderr, "\n");
        va_end(args);
        abort();
}

int x, y;

int width, height;
png_byte color_type;
png_byte bit_depth;

png_structp png_ptr;
png_infop info_ptr;
int number_of_passes;
png_bytep * row_pointers, *final_row_pointers;
int *foreground;
int my_cpu_id,numthreads;

int *final_tag, *final_level, tagindex ;
int *levelg ;
graph *G;

void alarm_wakeup (int j)
{
	int i,v;
	char *filename = "dummy.png" ;
	signal(SIGALRM,alarm_wakeup);
	printf("%d sec up partner, Wakeup!!!\n",INTERVAL);

		//setting the color
            for(i=0;i<G->nv;i++){
                if((final_level[i] == -1) && (levelg[i]) != -1){
                        final_level[i] = levelg[i];
                        final_tag[i] = tagindex;
                }
            }
            


		for (i=0; i<height; i++) {
			png_byte* row = final_row_pointers[i];
			for (v=0; v<width; v++) {
				//if(final_tag[i*width+v] != 0)
				//printf("final tag of %d is %d",i*width+v, final_tag[i*width+v]);
				row[v] = final_tag[i*width+v];
				//row[v] = 100;
			}
		}

		write_png_file(filename);
	
   exit(0);
}

void read_png_file(char* file_name)
{
        char header[8];    // 8 is the maximum size that can be checked

        /* open file and test for it being a png */
        FILE *fp = fopen(file_name, "rb");
        if (!fp)
                abort_("[read_png_file] File %s could not be opened for reading", file_name);
        fread(header, 1, 8, fp);
        if (png_sig_cmp(header, 0, 8))
                abort_("[read_png_file] File %s is not recognized as a PNG file", file_name);


        /* initialize stuff */
        png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);

        if (!png_ptr)
                abort_("[read_png_file] png_create_read_struct failed");

        info_ptr = png_create_info_struct(png_ptr);
        if (!info_ptr)
                abort_("[read_png_file] png_create_info_struct failed");

        if (setjmp(png_jmpbuf(png_ptr)))
                abort_("[read_png_file] Error during init_io");

        png_init_io(png_ptr, fp);
        png_set_sig_bytes(png_ptr, 8);

        png_read_info(png_ptr, info_ptr);

        width = png_get_image_width(png_ptr, info_ptr);
        height = png_get_image_height(png_ptr, info_ptr);


        color_type = png_get_color_type(png_ptr, info_ptr);
        bit_depth = png_get_bit_depth(png_ptr, info_ptr);

        number_of_passes = png_set_interlace_handling(png_ptr);
        png_read_update_info(png_ptr, info_ptr);


        /* read file */
        if (setjmp(png_jmpbuf(png_ptr)))
                abort_("[read_png_file] Error during read_image");

        row_pointers = (png_bytep*) malloc(sizeof(png_bytep) * height);
        for (y=0; y<height; y++)
                row_pointers[y] = (png_byte*) malloc(png_get_rowbytes(png_ptr,info_ptr));

	final_row_pointers = (png_bytep*) malloc(sizeof(png_bytep) * height);
        for (y=0; y<height; y++)
                final_row_pointers[y] = (png_byte*) malloc(png_get_rowbytes(png_ptr,info_ptr));

        png_read_image(png_ptr, row_pointers);

        fclose(fp);
}


void write_png_file(char* file_name)
{
        /* create file */
        FILE *fp = fopen(file_name, "wb");
        if (!fp)
                abort_("[write_png_file] File %s could not be opened for writing", file_name);


        /* initialize stuff */
        png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);

        if (!png_ptr)
                abort_("[write_png_file] png_create_write_struct failed");

        info_ptr = png_create_info_struct(png_ptr);
        if (!info_ptr)
                abort_("[write_png_file] png_create_info_struct failed");

        if (setjmp(png_jmpbuf(png_ptr)))
                abort_("[write_png_file] Error during init_io");

        png_init_io(png_ptr, fp);


        /* write header */
        if (setjmp(png_jmpbuf(png_ptr)))
                abort_("[write_png_file] Error during writing header");

        png_set_IHDR(png_ptr, info_ptr, width, height,
                     bit_depth, color_type, PNG_INTERLACE_NONE,
                     PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);

        png_write_info(png_ptr, info_ptr);


        /* write bytes */
        if (setjmp(png_jmpbuf(png_ptr)))
                abort_("[write_png_file] Error during writing bytes");

        png_write_image(png_ptr, final_row_pointers);


        /* end write */
        if (setjmp(png_jmpbuf(png_ptr)))
                abort_("[write_png_file] Error during end of write");

        png_write_end(png_ptr, NULL);

        /* cleanup heap allocation */
        for (y=0; y<height; y++){
                free(row_pointers[y]);
		free(final_row_pointers[y]);
         }
        free(row_pointers);
        free(final_row_pointers);

        fclose(fp);
}


int create_png_graph(int **tailp, int **headp)
{
        /*if (png_get_color_type(png_ptr, info_ptr) == PNG_COLOR_TYPE_RGB)
                abort_("[process_file] input file is PNG_COLOR_TYPE_RGB but must be PNG_COLOR_TYPE_RGBA "
                       "(lacks the alpha channel)");

        if (png_get_color_type(png_ptr, info_ptr) != PNG_COLOR_TYPE_RGBA)
                abort_("[process_file] color_type of input file must be PNG_COLOR_TYPE_RGBA (%d) (is %d)",
                       PNG_COLOR_TYPE_RGBA, png_get_color_type(png_ptr, info_ptr));*/

	int thisnode, nedges = 0;
	int max_edges = 100000000;
	*tailp = (int *) calloc(max_edges, sizeof(int));
    	*headp = (int *) calloc(max_edges, sizeof(int));
	//FILE *fp;
	
	foreground = (int *) calloc(width*height,sizeof(int));

	//fp = fopen ("png.temp","w");
	// y == row
	// x == column
        for (y=0; y<height; y++) {
                png_byte* row = row_pointers[y];
                for (x=0; x<width; x++) {
                        //png_byte* ptr = &(row[x*4]);
			png_byte* ptr = &(row[x]);
                        //printf("Pixel at position [ %d - %d ] has values: %d \n",x, y, ptr[0]);

                        /* set red value to 0 and green value to the blue one */
			if (nedges >= max_edges) {
			    printf("Limit of %d edges exceeded.\n",max_edges);
			    exit(1);
			}
			if(ptr[0]==255){
				thisnode = y*width+x;
				foreground[thisnode] = 1;
				if((y==0)&&(x==0)){
					png_byte* row1 = row_pointers[y+1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* ptr2 = &(row[x+1]);
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+1;}
				}
				else if((y==0)&&(x!=0)&&(x!=width-1)){
					png_byte* row1 = row_pointers[y+1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* ptr2 = &(row[x+1]);
					png_byte* ptr3 = &(row[x-1]);
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+1;}
					if(ptr3[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-1;}

				}
				else if((y==0)&&(x==width-1)){
					png_byte* row1 = row_pointers[y+1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* ptr2 = &(row[x-1]);
					
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-1;}

				}
				else if((y!=0)&&(x==0)&&(y!=height-1)){
					png_byte* row1 = row_pointers[y-1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* row2 = row_pointers[y+1];
					png_byte* ptr2 = &(row2[x]);
					png_byte* ptr3 = &(row[x+1]);
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+width;}
					if(ptr3[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+1;}

				}
				else if((y!=0)&&(x==width-1)&&(y!=height-1)){
					png_byte* row1 = row_pointers[y-1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* row2 = row_pointers[y+1];
					png_byte* ptr2 = &(row2[x]);
					png_byte* ptr3 = &(row[x-1]);
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+width;}
					if(ptr3[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-1;}

				}
				else if((y==height-1)&&(x==0)){
					png_byte* row1 = row_pointers[y-1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* ptr2 = &(row[x+1]);
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+1;}
				}
				else if((y==height-1)&&(x!=0)&&(x!=width-1)){
					png_byte* row1 = row_pointers[y-1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* ptr2 = &(row[x-1]);
					png_byte* ptr3 = &(row[x+1]);
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-1;}
					if(ptr3[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+1;}
				}
				else if((y==height-1)&&(x==width-1)){
					png_byte* row1 = row_pointers[y-1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* ptr2 = &(row[x-1]);
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-1;}
				}
				else{
					png_byte* row1 = row_pointers[y-1];
					png_byte* ptr1 = &(row1[x]);
					png_byte* row2 = row_pointers[y+1];
					png_byte* ptr2 = &(row2[x]);
					png_byte* ptr3 = &(row[x-1]);
					png_byte* ptr4 = &(row[x+1]);
					if(ptr1[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-width;}
					if(ptr2[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+width);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+width;}
					if(ptr3[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode-1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode-1;}
					if(ptr4[0]==255){
						//fprintf(fp,"%d %d\n",thisnode,thisnode+1);
						(*tailp)[nedges] = thisnode;
						(*headp)[nedges++] = thisnode+1;}
				}
			}
                        	
                        //ptr[1] = ptr[2];
                }
        }
	//fclose(fp);
	return nedges;
}


int read_edge_list (int **tailp, int **headp) {
    int max_edges = 100000000;
    int nedges, nr, t, h;
    FILE *fp;
    fp = fopen("png.temp","r");
    *tailp = (int *) calloc(max_edges, sizeof(int));
    *headp = (int *) calloc(max_edges, sizeof(int));
    nedges = 0;
    nr = fscanf(fp,"%i %i",&t,&h);
    while (nr == 2) {
	if (nedges >= max_edges) {
	    printf("Limit of %d edges exceeded.\n",max_edges);
	    exit(1);
	}
	(*tailp)[nedges] = t;
	(*headp)[nedges++] = h;
	nr = fscanf(fp,"%i %i",&t,&h);
    }
    fclose(fp);
    return nedges;
}


graph * graph_from_edge_list (int *tail, int* head, int nedges) {
    graph *G;
    int i, e, v, maxv;
    G = (graph *) calloc(1, sizeof(graph));
    G->ne = nedges;
    maxv = 0;

    // count vertices
    for (e = 0; e < G->ne; e++) {
	if (tail[e] > maxv) maxv = tail[e];
	if (head[e] > maxv) maxv = head[e];
    }
    G->nv = maxv+1;
    G->nbr = (int *) calloc(G->ne, sizeof(int));
    G->firstnbr = (int *) calloc(G->nv+1, sizeof(int));

    // count neighbors of vertex v in firstnbr[v+1],
    for (e = 0; e < G->ne; e++) G->firstnbr[tail[e]+1]++;

    // cumulative sum of neighbors gives firstnbr[] values
    for (v = 0; v < G->nv; v++) G->firstnbr[v+1] += G->firstnbr[v];

    // pass through edges, slotting each one into the CSR structure
    for (e = 0; e < G->ne; e++) {
	i = G->firstnbr[tail[e]]++;
	G->nbr[i] = head[e];
    }
    // the loop above shifted firstnbr[] left; shift it back right
    for (v = G->nv; v > 0; v--) G->firstnbr[v] = G->firstnbr[v-1];
    G->firstnbr[0] = 0;
    return G;
}


void print_CSR_graph (graph *G) {
    int vlimit = 20;
    int elimit = 50;
    int e,v;
    printf("\nGraph has %d vertices and %d edges.\n",G->nv,G->ne);
    printf("firstnbr =");
    if (G->nv < vlimit) vlimit = G->nv;
    for (v = 0; v <= vlimit; v++) printf(" %d",G->firstnbr[v]);
    if (G->nv > vlimit) printf(" ...");
    printf("\n");
    printf("nbr =");
    if (G->ne < elimit) elimit = G->ne;
    for (e = 0; e < elimit; e++) printf(" %d",G->nbr[e]);
    if (G->ne > elimit) printf(" ...");
    printf("\n\n");
}

void bfs_cpu_hybrid (int s, graph *G, int **levelp, int *nlevelsp, 
	int **levelsizep, int **parentp) {
    int *level, *levelsize;
    int thislevel;
    int **nnodearray,*cnodearray,*nsize, lqueuesize=10000000;
    int i, v, w, e,j,k,index;
    int size=0;
    int state = 0;
    int T1, T2, T3,alpha,beta;
    level = *levelp = (int *) calloc(G->nv, sizeof(int));
    levelsize = *levelsizep = (int *) calloc(G->nv, sizeof(int));

    // initially, queue is empty, all levels and parents are -1
    T1 = 64;
    if((G->nv)*0.01 > 262144) {
	T2 = (G->nv)*0.01;
    } else {
	T2 = 262144;
    }
    T3 = 2048;
    alpha = 2;
    beta = 2;

    numthreads=omp_get_max_threads();
 
    cnodearray = (int *) malloc(lqueuesize*sizeof(int));
    nnodearray = (int **) malloc(numthreads*sizeof(int *));
    for (i = 0; i < numthreads; i++){
	nnodearray[i] = (int *) malloc(lqueuesize*sizeof(int));
    }
   
    nsize = (int *) calloc(numthreads,sizeof(int));
#pragma omp parallel for private(v)
	for (v = 0; v < G->nv; v++) {
	    level[v] = -1;
	}

    // assign the starting vertex level 0 and put it on the queue to explore
    thislevel = 0;
    level[s] = 0;
    size=0;
    levelsize[0] = 1;
    cnodearray[0] = s;

    // loop over levels, then over vertices at this level, then over neighbors
    //printf("T1: %d , T2: %d, T3: %d \n",T1, T2, T3);
    while (levelsize[thislevel] > 0) {
	size=0;
	levelsize[thislevel+1] = 0;
	switch(state) {
	    case 0:
		if (levelsize[thislevel] > T1)
		    state =1;
		break;
	    case 1:
		if (levelsize[thislevel] > T2 )
		    state =2;
		if (levelsize[thislevel] > alpha*levelsize[thislevel-1] && levelsize[thislevel] > T3 )
		    state =2;
		break;
	    case 2:
		state = 3;
		break;
	    case 3:
		if (levelsize[thislevel] < T2)
		    state =4;
		if (levelsize[thislevel] < beta*levelsize[thislevel-1])
		    state =4;
		break;
	    case 4:
		if(levelsize[thislevel] > T1)
		    state = 1;
		else
		    state = 0;
		break;
	}
	//printf("state: %d , levelnum: %d, levelsize: %d \n",state, thislevel, levelsize[thislevel]);

	if (state == 2 || state == 3 ) {
#pragma omp parallel for private(i,e,w) reduction(+:size)
	    for (i = 0; i < G->nv; i++) {
		if(level[i]== thislevel) {
		    for (e = G->firstnbr[i]; e < G->firstnbr[i+1]; e++) {
			w = G->nbr[e];          // w is the current neighbor of v
			if (level[w] == -1) {   // w has not already been reached
			    level[w] = thislevel+1;
			    size++;
			}
		    }
		}
	    }
	} 
	else if (state == 4) {
#pragma omp parallel private(i,e,w,my_cpu_id)
	    {
		my_cpu_id=omp_get_thread_num();
# pragma omp for schedule(dynamic,CHUNKSIZE) nowait
		for (i=0; i < G->nv; i++){
		    if(level[i]== thislevel) {
			for (e = G->firstnbr[i]; e < G->firstnbr[i+1]; e++) {
			    w = G->nbr[e];          // w is the current neighbor of v

			    if (level[w] == -1) {   // w has not already been reached
				level[w] = thislevel+1;
				nnodearray[my_cpu_id][nsize[my_cpu_id]] = w;
				nsize[my_cpu_id]++;		
			    }
			}
		    }
		}            
	    }
	    size = 0;
	    index=0;
	    for (i = 0; i<numthreads ; i++){
		size = size + nsize[i];
		for(j=0;j<nsize[i];j++){
		    cnodearray[index] = nnodearray[i][j];
		    index++;
		}
		nsize[i] = 0;
	    }
	}
	else if (state == 1) {
#pragma omp parallel private(i,e,w,v,my_cpu_id)
	    {
		my_cpu_id=omp_get_thread_num();
# pragma omp for schedule(dynamic,CHUNKSIZE) nowait
		for (i=0; i<levelsize[thislevel]; i++){
		    v = cnodearray[i];
		    //printf("v=%d\n",v);
		    for (e = G->firstnbr[v]; e < G->firstnbr[v+1]; e++) {
			w = G->nbr[e];          // w is the current neighbor of v

			if (level[w] == -1) {   // w has not already been reached
			    level[w] = thislevel+1;
			    nnodearray[my_cpu_id][nsize[my_cpu_id]] = w;
			    nsize[my_cpu_id]++;		
			}
		    }
		}            
	    }
	    size = 0;
	    index=0;
	    for (i = 0; i<numthreads ; i++){
		size = size + nsize[i];
		for(j=0;j<nsize[i];j++){
		    cnodearray[index] = nnodearray[i][j];
		    index++;
		}
		nsize[i] = 0;
	    }
	}
	else {
	    for (i=0; i<levelsize[thislevel]; i++){
		v = cnodearray[i];
		for (e = G->firstnbr[v]; e < G->firstnbr[v+1]; e++) {
		    w = G->nbr[e];          // w is the current neighbor of v
		    if (level[w] == -1) {   // w has not already been reached
			level[w] = thislevel+1;
			nnodearray[0][nsize[0]] = w;
			nsize[0]++;		
		    }
		}
	    }            
	    size = 0;
	    index=0;
	    size = size + nsize[0];
	    for(j=0;j<nsize[0];j++){
		cnodearray[index] = nnodearray[0][j];
		index++;
	    }
	}
	thislevel = thislevel+1;
	levelsize[thislevel]=size;
    }
    *nlevelsp = thislevel;
    //free(queue);
}


void bfs_read (int s, graph *G, int **levelp, int *nlevelsp, 
	int **levelsizep, int **parentp) {
    int *level, *levelsize, *parent;
    int thislevel;
    int *queue, back, front;
    int i, v, w, e;
    int chunk,size;
    int extra_part;
    level = *levelp = (int *) calloc(G->nv, sizeof(int));
    levelsize = *levelsizep = (int *) calloc(G->nv, sizeof(int));
    //parent = *parentp = (int *) calloc(G->nv, sizeof(int));
    //queue = (int *) calloc(G->nv, sizeof(int));

    // initially, queue is empty, all levels and parents are -1
    //back = 0;   // position next element will be added to queue
    //front = 0;  // position next element will be removed from queue
//#pragma omp parallel private(v)
//    {
//# pragma omp for nowait
#pragma omp parallel for private(v)
	for (v = 0; v < G->nv; v++) {
	    //parent[v] = -1;
	    level[v] = -1;
	}
//    }

    // assign the starting vertex level 0 and put it on the queue to explore
    thislevel = 0;
    level[s] = 0;
    size=0;
    levelsize[0] = 1;
  // loop over levels, then over vertices at this level, then over neighbors
    while (levelsize[thislevel] > 0) {
        size=0;
#pragma omp parallel for private(i,e,w) reduction(+:size)
        for (i = 0; i < G->nv; i++) {
            if(level[i]== thislevel) {
                for (e = G->firstnbr[i]; e < G->firstnbr[i+1]; e++) {
                    w = G->nbr[e];          // w is the current neighbor of v
                    if (level[w] == -1) {   // w has not already been reached
                        level[w] = thislevel+1;
                        size++;
                    }
                }
            }
        }
        thislevel = thislevel+1;
        levelsize[thislevel]=size;
    }
    *nlevelsp = thislevel;

}


void bfs_queue (int s, graph *G, int **levelp, int *nlevelsp, 
	int **levelsizep, int **parentp) {

    int *level, *levelsize, *parent;
    int thislevel;
    //int *queue, back, front,*lqueue, lback, lfront, lqueuesize=1000000, levelval;
    int **nnodearray,*cnodearray,*nsize, lqueuesize=10000000;
    int i, v, w, e,j,k,index;
    int size = 0;;
    level = *levelp = (int *) calloc(G->nv, sizeof(int));
    levelsize = *levelsizep = (int *) calloc(G->nv, sizeof(int));
    
    numthreads=omp_get_max_threads();
    //printf("numthreads in bfs = %d\n",numthreads);
    
    cnodearray = (int *) malloc(lqueuesize*sizeof(int));
    nnodearray = (int **) malloc(numthreads*sizeof(int *));
    for (i = 0; i < numthreads; i++){
	nnodearray[i] = (int *) malloc(lqueuesize*sizeof(int));
        //nnodearray[i] = (int *) malloc(lqueuesize*sizeof(int));
    }
    
    
    nsize = (int *) calloc(numthreads,sizeof(int));
    //nsize = (int *) calloc(numthreads,sizeof(int));

    //#pragma omp parallel private(v)
    {
    //# pragma omp for nowait
    for (v = 0; v < G->nv; v++) {
        //parent[v] = -1;
	level[v] = -1;
    }
    }

    // assign the starting vertex level 0 and put it on the queue to explore
    thislevel = 0;
    level[s] = 0;
    size=0;
    levelsize[0] = 1;
    cnodearray[0] = s;

	
    // loop over levels, then over vertices at this level, then over neighbors
    while (levelsize[thislevel] > 0) {
	levelsize[thislevel+1] = 0;
	#pragma omp parallel private(i,e,w,v,my_cpu_id)
	{
	    my_cpu_id=omp_get_thread_num();
	    # pragma omp for schedule(dynamic,CHUNKSIZE) nowait
	    for (i=0; i<levelsize[thislevel]; i++){
		v = cnodearray[i];
		//printf("v=%d\n",v);
		    for (e = G->firstnbr[v]; e < G->firstnbr[v+1]; e++) {
			w = G->nbr[e];          // w is the current neighbor of v
	
			if (level[w] == -1) {   // w has not already been reached
			    level[w] = thislevel+1;
			    //final_tag[w] = tagindex;
			    nnodearray[my_cpu_id][nsize[my_cpu_id]] = w;
			    nsize[my_cpu_id]++;		
			}
		    }
		//}
	    }            
	}
	//printf("this level = %d\n",thislevel);
	size = 0;
	index=0;
	for (i = 0; i<numthreads ; i++){
		size = size + nsize[i];
		for(j=0;j<nsize[i];j++){
			cnodearray[index] = nnodearray[i][j];
			index++;
		}
		nsize[i] = 0;
	}
	
	thislevel = thislevel+1;
	levelsize[thislevel]=size;
	//printf("this size = %d\n",size);
    }
	*nlevelsp = thislevel;
	//free(queue);
}

void bfs_ser (int s, graph *G, int **levelp, int *nlevelsp, 
	int **levelsizep, int **parentp) {
    int *level, *levelsize, *parent;
    int thislevel;
    int *queue, back, front;
    int i, v, w, e;
    level = *levelp = (int *) calloc(G->nv, sizeof(int));
    levelsize = *levelsizep = (int *) calloc(G->nv, sizeof(int));
    parent = *parentp = (int *) calloc(G->nv, sizeof(int));
    queue = (int *) calloc(G->nv, sizeof(int));

    // initially, queue is empty, all levels and parents are -1
    back = 0;   // position next element will be added to queue
    front = 0;  // position next element will be removed from queue
    for (v = 0; v < G->nv; v++) level[v] = -1;
    for (v = 0; v < G->nv; v++) parent[v] = -1;

    // assign the starting vertex level 0 and put it on the queue to explore
    thislevel = 0;
    level[s] = 0;
    levelsize[0] = 1;
    queue[back++] = s;

    // loop over levels, then over vertices at this level, then over neighbors
    while (levelsize[thislevel] > 0) {
	levelsize[thislevel+1] = 0;
	for (i = 0; i < levelsize[thislevel]; i++) {
	    v = queue[front++];       // v is the current vertex to explore from
	    for (e = G->firstnbr[v]; e < G->firstnbr[v+1]; e++) {
		w = G->nbr[e];          // w is the current neighbor of v
		if (level[w] == -1) {   // w has not already been reached
		    parent[w] = v;
		    level[w] = thislevel+1;
		    //printf("w = %d level = %d\n",w,level[w]);	
		    levelsize[thislevel+1]++;
		    queue[back++] = w;    // put w on queue to explore
		}
	    }
	}
	thislevel = thislevel+1;
    }
    *nlevelsp = thislevel;
    free(queue);
}


int main (int argc, char* argv[]) {
    //graph *G;
    int *levelsize, *parent , finished = 0;
    //int tagindex ;
    int *tail, *head;
    int nedges;
    int nlevels;
    int startvtx;
    int i, v, reached;
    double time_elapsed;
    double start,end;
    //int *final_level;
    //clock_t start=clock();

    struct itimerval tout_val;

    tout_val.it_interval.tv_sec = 0;
    tout_val.it_interval.tv_usec = 0;
    tout_val.it_value.tv_sec = INTERVAL; /* set timer for "INTERVAL (10) seconds */
    tout_val.it_value.tv_usec = US_INTERVAL;


    if (argc == 3) {
	startvtx = atoi (argv[1]);
    } else {
	printf("usage:   png_bfs_conn_comp <inputpng> <outputpng>\n");
	//printf("example: cat sample.txt | ./bfstest 1\n");
	exit(1);
    }


    omp_set_num_threads(16);
#ifdef _OPENMP
    //numthreads=omp_get_num_threads();
    numthreads=omp_get_max_threads();
    //numthreads = 1;
#else
    numthreads=1;
#endif
    //time_elapsed = omp_get_wtime();
    //start = time(NULL);

    // Read png file 
    read_png_file(argv[1]);


    nedges = create_png_graph(&tail, &head);

    //nedges = read_edge_list (&tail, &head);
    
    G = graph_from_edge_list (tail, head, nedges);
    free(tail);
    free(head);

    
    final_level = (int *) calloc(width*height, sizeof(int));
    final_tag = (int *) calloc(width*height, sizeof(int));
    tagindex = 1;//level

    for(i = 0; i<width*height;i++)
	final_level[i] = -1;

    print_CSR_graph (G);

    printf("numthreads = %d\n", numthreads);
    //printf("Starting vertex for BFS is %d.\n\n",startvtx);

    // Perform connected component analysis here

    for(i=0;i<G->nv;i++){
	if(foreground[i] == 1){
		final_level[i] = 0;
		final_tag[i] = 1;
		startvtx = i;
		break;
	}
    }

    
    setitimer(ITIMER_REAL, &tout_val,0);
    signal(SIGALRM,alarm_wakeup); /* set the Alarm signal capture */

    start=omp_get_wtime();
    
    //each loop for each component
    while(1){
	    finished = 1;
	    //final_level[startvtx] = 0;
	    //printf("Start vtx= %d\n",startvtx);
	    //bfs (startvtx, G, &level, &nlevels, &levelsize, &parent, final_level, final_tag, tagindex);
	    
	    //bfs_ser (startvtx, G, &level, &nlevels, &levelsize, &parent);
	    bfs_read (startvtx, G, &levelg, &nlevels, &levelsize, &parent);
	    //bfs_queue (startvtx, G, &level, &nlevels, &levelsize, &parent);
	    //bfs_cpu_hybrid (startvtx, G, &level, &nlevels, &levelsize, &parent);

	    # pragma omp parallel private(i) 
	    {
	    # pragma omp for nowait

	    //setting the color
	   
	    for(i=0;i<G->nv;i++){
		if((final_level[i] == -1) && (levelg[i]) != -1){
			final_level[i] = levelg[i];
			final_tag[i] = tagindex;
		}
	    }
	    }

	    //find startvertex for next iteration
	    for(i=startvtx;i<G->nv;i++){
		if((final_level[i] == -1) && (foreground[i]) == 1){
			tagindex = (tagindex + 1) % 255 ;
			//tagindex++;
			startvtx = i;
			finished = 0;
			break;
		}
	    }
	    if(finished == 1)
		break;
    }
    end=omp_get_wtime();

    for (i=0; i<height; i++) {
        png_byte* row = final_row_pointers[i];
        for (v=0; v<width; v++) {
		row[v] = final_tag[i*width+v];
	}
    }
    /*reached = 0;
    for (i = 0; i < nlevels; i++) reached += levelsize[i];
    printf("Breadth-first search from vertex %d reached %d levels and %d vertices.\n",startvtx, nlevels, reached);
    for (i = 0; i < nlevels; i++) printf("level %d vertices: %d\n", i, levelsize[i]);
    
    if (G->nv < 20) {
        printf("\n  vertex parent  level\n");
        //for (v = 0; v < G->nv; v++) printf("%6d%7d%7d\n", v, parent[v], level[v]);
        printf("\n  vertex level\n");
        for (v = 0; v < G->nv; v++) printf("%6d%7d\n", v, levelg[v]);
    }*/

    write_png_file(argv[2]);
   // write_png_file("dummy.png");
    printf("Number of connected components = %d \n",tagindex);
    printf("\n");
    printf("Time to run = %f.\n",end-start);
    //printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC);*/
}
